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Abstract. An abstract should be given 

his work reports numerical experiments intended 
to clarify the internal equilibration process in large 
molecules, following vibrational excitation. A model of 
amorphous and oxygenated hydrocarbon macromolecule 
(^500 atoms)-simulating interstellar dust- is built up by 
means of a chemical simulation code. Its structure is op- 
timized, and its normal modes determined. About 4.5 eV 
' of potential energy is then deposited locally by perturb- 
ing one of the C-H peripheral bonds, thus simulating the 
capture of a free H atom by a dangling C bond. The en- 
suing relaxation of the system is followed for up to 300 
ps, using a molecular mechanics code. When steady state 
is reached, spectra and time correlation functions of ki- 
netic energy and bond length fluctuations indicate that 
most normal modes have been activated, but the motion 
remains quasi-periodic or regular. By contrast, when the 
molecule is violently excited or embedded in a thermal 
bath (modelled by Langevin dynamics), the same markers 
clearly depict chaotic motions. Thus it appears that even 
such a large system of oscillators is unable to provide the 
equivalent of a thermal bath to any one of these, barring 
strong resonances between some of them. This conclusion 
is of consequence for the interpretation of astronomical 
infrared spectra. 

Collateral numerical experiments show that a) relax- 
ation times increase as perturbation energy decreases by 
spreading through the system; b) energy deposited in the 
highest- frequency modes does not relax preferentially into 
the lowest frequency-modes but follows specific paths de- 
termined by near resonances and coupling strength; c) en- 
ergy deposited in the lowest-frequency modes is able to 
flow up the whole frequency ladder, albeit less easily than 
in the opposite direction. 
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1. Introduction 

The present research was triggered by the question 
whether infrared chemiluminescencc was possible with 
large molecules. This question, in turn, arose from the 
need to interpret the flurry of accurate IR (infrared) spec- 
tra delivered by dedicated astrophysics satellites in the 
past two decades. It is indeed established that the ISM 
(InterStellar Medium) is interspersed with large molecules 
or small grains (dust), with a wide size distribution (up to 
1000 nm), and whose composition is essentially limited 
to hydrocarbons and silicates. We are interested, here, 
in that major family of amorphous hydrocarbons which 
lie so far from any illuminating star that they can only 
reach a low steady state temperature (<200 K) when the 
UV radiative power they absorb equals the IR radiative 
power they emit. Under such circumstances, the relatively 
high IR spectral intensities observed are incompatible with 
Planck's law at that temperature. A possible alternative 
excitation mechanism is chemiluminescence, based on the 
fact that, in such environments, the dominant material, 
hydrogen, remains in atomic form, which notably favours 
interaction with hydrocarbon dust. Such interaction is be- 
lieved to be the main cause of formation of molecular hy- 
drogen, which is the main constituent of the "clouds" that 
are detected within the ISM (Duley and Williams 1 ). When 
a free H atom hits an H atom that is attached to a C atom 
at the periphery of a dust particle, they may form such a 
strong bond that the nascent H2 molecule can fly away, 
leaving a dangling bond behind. Upon a later encounter 
with an H atom, the latter can be captured by the dangling 
bond to form a long-lived C-H bond. In such an event, an 
energy of order 4.5 eV (the bond dissociation energy) is 
deposited locally in the molecule. This work addresses the 
fate of that energy and of the resulting atomic vibrations: 
what vibrational modes are excited and can they survive 
energy redistribution long enough for them to finally relax 
through the relatively slow process of IR emission? 
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Of course, IR chcmiluminesccncc is known to occur 
upon reactions between small molecules (e.g. see Lee and 
Pimentel 2 ). Is this possible with larger molecules as well? 
and how much larger? More precisely, does statistical equi- 
librium set in, with equal probability for all states of equal 
energy (also referred to as complete randomization)? If it 
does, a temperature can be defined, which scales like the 
inverse of the number of atoms; then, for a dust particle 
larger than, say, 100 atoms, the temperature would be so 
low as to preclude detectable emission of all but the very 
lowest vibrational frequencies, contrary to observations. 
If, on the contrary, complete randomization does not ob- 
tain, then the motions of the various atoms retain at least 
some degree of coherence and the concept of temperature 
becomes irrelevant. The distribution of energy among the 
various modes is then totally different, being determined 
by classical ( or quantum, for more accuracy) mechanics: 
even the highest vibrational frequencies can be excited to 
a notable extent and IR chcmiluminescence becomes de- 
tectable. 

The RRKM theory of unimolecular reactions (see 
Gilbert and Smith 3 ) indeed assumes thermodynamical 
equilibrium, and is vindicated by innumerable experimen- 
tal data. However, these data are mostly about yields 
of particular unimolar reactions (dissociations) averaged 
over assemblies of molecules under finite pressure, so that 
vibrational excitation is effected collisionally and ran- 
domly, according to equilibrium thermodynamics. This, 
according to Freed 4 , is enough to ensure the validity of 
all the mathematical predictions of the RRKM theory, 
even assuming no randomization at all occurs within the 
molecule! Besides, more recent experiments have shown 
deviation from statistical behaviour, especially when the 
energy content is low (see McDonald and Marcus 5 ). Other 
types of experiments, especially those designed in pursuit 
of "selective chemistry", suggest that, under suitable con- 
ditions, the molecule can be prepared in an excited and 
coherent state, which can survive for quite long times (see 
Zewail 6 ). More recently, Williams and Leone followed the 
IR- radiative relaxation of gaseous naphtalene at very low 
pressures, after laser excitation of a C-H stretch vibra- 
tion; they showed that their spectra are not compatible 
with calculations based on thermodynamic equilibrium. 
Also, Deak et al. 8 excited the C-H stretch of nitromcthane 
molecules in a tiny segment of a jet, and observed the 
prompt (<2 ps) establishment of a nearly uniform dis- 
tribution among other vibrations, which prevailed for a 
time about equal to the collisional relaxation time of the 
bath (30 to 50 ps). Afterwards, a new, Boltzmannian, 
distribution clearly set in, again suggesting that colli- 
sions where responsible for this to occur. Similar con- 
clusions were drawn by Stromberg et al. 9 , working on 
W(CO)6- Time and frequency resolved experiments reveal 
in fact three cases of IVR: dissipative (energy sharing), re- 
stricted, and no IVR at all (see, for instance, Zewail 10 ). 
Correspondingly, molecular motions range from statisti- 
cal, ergodic, chaotic, to quasi-periodic and periodic, or 
regular. 



Turning to theory, it is known that the normal modes 
of vibration of a system are orthogonal and, therefore, 
in principle do not interact. However, they do so in real 
life, because anharmonicity always provides some coupling 
between them, albeit weak. This coupling increases in 
strength as the molecule becomes larger, and the distance 
between modes correspondingly decreases in frequency 
space. There is therefore no doubt that any local distur- 
bance is bound to decay in some way, except, in principle, 
in the ideal limit of the weakest excitation of the purest 
normal mode. This type of reasoning led Fermi, in the 
30's, to conjecture that a statistical (Boltzmannian) equi- 
librium must ultimately set in (see Ford 11 ). He therefore 
initiated numerical simulations of the dynamics of long 
chains of oscillators, the outcome of which did not, how- 
ever, confirm his conjecture: hence the Fermi-Ulam-Pasta 
paradox (Fermi 12 ). Long afterwards it became clear that 
statistical intramolecular energy transfer is not ubiquitous 
(see Brumer 13 ) even with many degrees of freedom; con- 
versely, under suitable circumstances, "chaos" can occur 
even with a few degrees of freedom. Thus, the size of the 
system and the corresponding number of normal modes 
are not good indicators of fast randomization 14 . 

Finally, a hierarchy of statisticality (stochasticity) and 
randomness was established, with ergodic, mixing, K and 
Bcrnouilli systems in order of increasing randomness 15 . 
Conditions favourable to chaos and randomness were rec- 
ognized (see Ford 11 ' 16 ): existence and overlapping of non- 
linear resonances (Sn^^e), high energy content, brief 
random excitations (as in a perfect gas or with billiard 
balls), etc. 

In terms of phase space, Boltzmann's ergodic hypoth- 
esis implies that all available phase space is uniformly 
visited by the system's trajectory during the observation 
time. This leads to the usual canonical prescription for de- 
riving time averages of physical quantities. It is found that 
energy is distributed among the degrees of freedom (here, 
the vibrational modes) according to Boltzmann's law with 
a single temperature common to all. While this delivers 
the least biased estimates of observables in the absence of 
enough information on the system, very many well docu- 
mented physical systems are not properly described in this 
way 17 . In these cases, the phase space is full of structures 
which induce "intramolecular bottlenecks" to the phase 
flow (see Uzer 18 for a review). Such systems are found to 
obey a restricted form of crgodicity, broken ergodicity , 
meaning that, in these cases, statisticality (randomness) 
obtains only in limited regions of phase space. 

For all these developments, theory can hardly predict 
which type of internal motion will prevail in systems with 
more than a few degrees of freedom. Surely, experiment 
remains the ultimate recourse. However, numerical mod- 
els of coupled oscillators have been of great help in high- 
lighting the role of different parameters (particularly ex- 
citation strength and bond anharmonicity), in assessing 
the degree of randomness, in interpreting experimental 
outcomes and in evaluating relaxation times (see Noid et 
al. 19 and Thompson and co-workers 20 ). With the advent of 
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highly performing chemical software, numerical modeling 
can now more closely simulate the motion of large, real- 
life molecules. This is the course followed in the present 
work to address the problem mentioned above. 

The main result of this research is as follows: in the 
type of structure considered here, and following the re- 
combination of an impinging H atom with a dangling C 
bond, the energy deposited is ultimately, and on average, 
roughly equipartitioned among all the vibrational modes. 
However, chaos does not ensue and the molecule remains 
in a nearly coherent state, i.e. a coherent superposition of 
normal vibration modes . A Fourier analysis of the time 
variations of the total electric dipole moment then pro- 
vides the spectral distribution of energy, which shows that 
chemiluminescence is indeed detectable. 

The originality of this work does not lie in the method 
used nor in the conclusion that statistical equilibrium does 
not always obtain after relaxation. Rather, it is in the 
combination of very large molecules (~ 500 atoms), amor- 
phous structure, chemical excitation, no collision at all 
before complete internal relaxation and very long compu- 
tation runs (up to 500 ps). Such a combination, suggested 
by astronomical conditions, is difficult to achieve in the 
laboratory and, to my knowledge, was never considered 
before. Moreover, the outcome of this study may seem 
counterintuitive to many who are not familiar with the 
FUP paradox, and therefore may fuel interesting debate. 

2. The chemical software 

The calculations were made using the Hypcrchem package 
released by Hypercube, Inc. in 1999, and implemented on a 
desk-top PC equipped with a Pentium (r)II processor (450 
MHz) and MMX (TM) technology. First, the molecule is 
built on the screen by picking atom after atom in a table 
of elements. The CC bondings are specified graphically: -, 
= , =. A computational code and options are then selected. 
The code is then asked to optimize the molecular geom- 
etry, i.e. determine the architecture which minimizes the 
total potential energy, Eo. A perturbation is then applied 
by displacing one or more atoms, or imparting a velocity to 
them, or creating a new bond (like in chemisorption) . This 
brings it up to an energy state E*. The subsequent mo- 
tion is then followed using a Molecular Dynamics routine 
based on the Born-Oppenhcimcr approximation: at each 
calculation step, the nuclei arc assumed to be clamped 
at their initial positions and the code determines the elec- 
tronic charge distribution, either by solving Schroedinger's 
equation ((a) below) or empirically ((b) below). Hence the 
forces acting on the nuclei and the new velocities of the 
latter, which determine the displacement for the next step. 
The step length is set at 10~ 16 s, much less than the short- 
est characteristic period of the system. To save time, the 
electric dipole moment is computed only once every cy- 
cle of 30 dynamics steps, still ensuring correct frequency 
analysis up to 3000 cm -1 . 

For a frequency resolution of 30 cm -1 , a Molecular 
Dynamics run must span 1 ps, i.e. 10000 steps or 300 cy- 



cles. However, much longer times are needed for transients 
to decay. The coordinates and velocity of each nucleus, and 
hence the potential and kinetic energies are recomputed 
at every step. The total energy is conserved all along, as 
it should be, even for the longest calculation. The energy 
content, or perturbation energy is Ep= E*-Eo . The accu- 
racy and speed of this computation depend on the choice 
of the code; while most state-of-the-art chemical codes are 
available with this package, only two of them were used in 
the present work: 

a) the semi-empirical PM3 code, developed by J. 
Stewart 21 along the same line as AMI. The semi-empirical 
methods use a rigorous quantum-mechanical formalism 
combined with empirical parameters obtained from com- 
parison with experimental results. PM3 incorporates a 
much larger number and wider variety of experimental 
data than AMI. These codes compute approximate so- 
lutions of Schroedinger's equation, and use experimental 
data only when the quantum mechanical calculations are 
too difficult or too lengthy. This sometimes makes them 
more accurate than poor ab initio methods, and they are 
always faster. PM3 also provides for a direct calculation 
of the normal mode spectrum, which does not require 
running the Molecular Dynamics. This is based on the 
quadratic term of the interatomic potentials and only in- 
volves matrix manipulation. The result applies, more gen- 
erally, to weak excursions from the equilibrium geometry. 
The spectrum is discrete, i.e. composed of separate vertical 
bars of length equal to the "IR intensity" , A, of the struc- 
ture at the corresponding wavenumber. In real life, this 
must be corrected for broadening effects: thermal fluctu- 
ations, rotations, and slight differences between particles 
of the same family. For large particles, the first two effects 
are negligible in the astrophysical context. 

The accuracy of the mode wavenumbers has been dis- 
cussed at length 22 : although it depends on the type of 
vibration and on the structure, it is considered to be of 
order 10 cm" 1 for the C-H stretch and 100-200 for the rest 
of the spectrum. 

An interesting feature of this code is the possibility of 
animating the atoms on the screen in any selected nor- 
mal mode. This is very useful in the process of tailoring 
the composition and structure of the particle to fit the 
astronomical spectra as best as possible. 

This code is very much faster than ab initio codes 
(which are out of the question here). However, for a 
molecule of 100 atoms, each Molecular Dynamics step re- 
quires ~30 s computation time, which amounts to about 
100 cycles a day or 1 month for a run of minimum length. 
This code was used in the present work to determine, by 
cut-and-try, a structure capable of displaying the main IR 
features observed in the sky; this structure is shown in fig. 
1 (Papoular 23 , 24 ). Note that it is 3-D, heterogeneous, dis- 
ordered, amorphous and closely resembles models of coal 25 
and kerogen 26 , a material also found in meteorites. In this 
study, we have used larger molecules built up by assem- 
bling randomly part or totality of fig. 1 (101 atoms) so as 
to retain the same chemical characters. Most of the results 
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below refer to the largest of these molecules, designated 
by LM and consisting of 493 atoms. This is too large to 
be conveniently handled with semi-empirical codes. 

b) The Molecular Mechanics code MM+. This is based 
on the Allinger 27 MM2 code. Molecular Mechanics codes, 
like MM2, Amber, Charmm,... approximate the required 
generalized potential energy surface (force fields) by an- 
alytical functions which are deduced from, and subse- 
quently parametrized via, accurate calculations and ex- 
perimental results. The force field contains terms associ- 
ated each with bond stretching, angle (in-plane) bending, 
torsional (out-of-plane) bending, stretch-bend cross terms, 
twisting, van der Waals forces and other terms of chemical 
significance. Anharmonicity is arbitrarily built in various 
terms; it is cubic for bond stretches. 

MM+ uses the latest MM2 parameters and atom types 
(designed primarily for small organic molecules), comple- 
mented, where necessary, by specific calculations. When 
the "atomic charge" option is selected, the charge on ev- 
ery atom is invariable and determined empirically by the 
type of atom and bonding. As a result, the vibration fre- 
quencies derived from MM+ calculations may be in error 
by up to 200 cm -1 . Moreover, the absolute value of the 
electric dipole moment is not reliable. However, this code 
is very much faster than the semi-empirical codes: about 
20000 cycles a day for a particle of 100 atoms. It is, there- 
fore, highly useful in clarifying the general trends of the 
spectra as a function of parameters such as time, energy 
content, structure, etc. Such trends, which are our main 
interest here, are not likely to depend upon the exact nu- 
merical values of vibration frequencies in large molecules 
like our's. This code was therefore used extensively for 
our present purposes; since it is not suited to simulating 
the hydrogen recombination with the molecule, the lat- 
ter is approximately simulated by starting with the model 
particle at rest in its ground state, Eo, as above . Then, 
instead of breaking a C-H bond, the bond is simply altered 
in length and direction on the screen, by cut and try, until 
the new potential energy of the particle is E*~Eo+4.5 eV. 
The Molecular Dynamics is launched thereafter. Under the 
perturbation just described, the real molecule is expected 
to remain in its electronic ground state. Energy is con- 
served to at least 3 significant digits and much better on 
average over a few steps. 

The relevance and efficiency of force fields and semi- 
empirical methods are now well established in spite of their 
neglect of quantum effects. The rationale behind their ex- 
tensive use lies in the Ehrenfest theorem (see 34 for a dis- 
cussion). 

3. The approach to steady state 

The quantitative criterion adopted here for stationarity is 

R=M°- 5 ( r/(E fcm } (1) 

where M = 3N — 6 is the number of degrees of freedom 
for TV atoms, (E^n) is the average total kinetic energy, and 




Fig. 1. The optimized structure of the building block (101 
atoms) of the reference macromolecule, LM (493 atoms). 
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Fig. 2. The reduced standard deviation, R=M 5 a/ (E kin } 
, of the total kinetic energy content, Efc ira , of the molecule, 
as a function of time after local injection of 4.5 eV of 
chemical energy; a is the std. dev. and M the number of 
normal modes. The dashed curve is a rough fit with the 
extended exponential, l+30/exp(t/0.8) ' 4 . 



er, its standard deviation. As a function of time, this tends 
towards a horizontal asymptote. This is shown in fig. 2 
for our largest molecule (LM, 493 atoms). The relaxation 
is initially very fast (r~0.5 ps) but later looks more like 
an extended exponential (dashed curve). Ultimately, R«l 
and all the degrees of freedom appear to have been acti- 
vated. 

The code allows us to follow the "diffusion" of en- 
ergy in time and space, from the perturbed bond into the 
molecule ; this is done by defining a circle/sphere centered 
on the captured H atom, "selecting" the atoms within 
and asking for their total kinetic energy; this procedure 
is repeated for any desired radius and instant along the 
calculation. Concurrently, autocorrelation and FFT (Fast 
Fourier Transform) can be obtained from the time evolu- 
tion of the position/velocity of any atom and of the length 
or orientation of any bond, revealing the nature of the dy- 
namics. In this way, it is found that 
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Fig. 3. The relaxation of the cumulative spectral kinetic 
energy distribution, Y^kin{v) , where Efcj„(z/) is the 
Fourier transform of Efc, n (t); normalized to 1 at the high- 
est computed frequency. Curves 1 to 7 were recorded suc- 
cessively after delays of 0.41, 1.6, 5, 55, 71, 121, and 221 
ps, respectively, starting at the beginning of the run. 

1) the deposited chemical (or potential) energy is 
promptly transformed into stretch (^3000 cm _1 )and bend 
(~1500 cm -1 ) vibrations of the central C-H bond; these 
motions are initially highly non-linear and carry strong 
harmonics; 

2) spatially extended low-frequency vibrations quickly 
invade the whole molecule while high-frequency vibra- 
tions lag behind because they are more localized in space 
and preferentially exchange energy between neighbouring 
atoms; 

3) overall, low frequencies grow at the expense of high 
frequencies, until all modes are activated and rough energy 
cquipartition obtains between all atoms; 

4) steady state is essentially established at about 100 

ps; 

One would also wish to know the energy in each nor- 
mal mode. However, a normal mode analysis is hopeless 
because of the system size. A less rigorous, but simpler 
and quite helpful, approach consists in Fourier analyzing 
the total kinetic energy, Efcj„, over an arbitrary period 
(here from 0.4 to 16 ps). Because of spectral congestion, 
this " energy" spectrum is inconvenient for illustrative pur- 
poses. It is better to use the cumulative spectral kinetic 
energy distribution, X)Efcm(^), which is obtained by sum- 
ming up the spectrum from to v, the abscissa. The result 
is drawn in fig. 3 for several successive times during the 
transient; the curves are all normalized to 1 at the high- 
est frequency. Again because of spectral congestion, they 
seem to be roughly monotonous on this scale, except for 
the isolated C-H stretch band near 6000cm -1 . Since Efcm 
is a quadratic form in the coordinates or momenta, the 
frequencies are, in fact, double the vibration frequencies. 
In case of a purely harmonic motion, one also expects sum 
and difference frequencies. These are washed out, however, 



in the presence of strong dephasing between all frequency 
components, which seems to be the case here, making for 
an easier interpretation. 

Points 3) and 4) made above are synthetically illus- 
trated in fig. 3. It appears that the spectrum is essentially 
invariable after a few picoseconds, although some "redden- 
ing" goes on for a much longer time. It is remarkable, at 
this point, that the energy initially concentrated in the C- 
H stretch flows roughly from high to low frequencies: the 
main relaxation path is not direct to the "bath" frequen- 
cies. Neither is it a strict "vibrational cascade" down a lad- 
der of levels: significant amounts of intermediate and low 
frequency energy are observed at the earliest stage; more- 
over, the relaxation time roughly increases with decreasing 
mode frequency All this was already deduced from exper- 
iments on small molecules in liquid environment (see Deak 
et al. 1999 8 , Sibert and Rey 8 ). In the present study, the 
system is isolated so the energy decay of the C-H stretch 
modes (5000 to 6000 cm -1 in fig. 3) is only partial; as ex- 
pected from anharmonic oscillators, it is accompanied by 
line narrowing and blue shift . 

4. The steady state 

We now take up the detail character of the motion in 
steady state. There are a number of different features 
which distinguish the various types of molecular motions 
(Noid et al. 19 ). Most are associated with phase space coor- 
dinates; these are well suited to theoretical studies but can 
hardly be accessed in the present case. More manageable 
are the methods which involve looking at the frequency 
spectra or correlation functions of dynamical variables. 
As noted by Noid et al., the frequency spectrum of a dy- 
namical variable, for quasi-periodic (regular) motion, con- 
sists of a series of sharp lines, while the chaotic frequency 
spectrum has broad lines forming bands near the location 
of the sharp lines of the corresponding regular spectrum. 
The transition from regular to chaotic is also characterized 
by the continuum between bands becoming less and less 
smooth and more and more rugged (grassy). In principle, 
for a given perturbation, every normal mode shows up, to 
a smaller or larger extent, in the fluctuations of position 
and velocity of each atom, as well as of the length of each 
bond, depending on the type of atom or bond and immedi- 
ate environment. C-C bonds were found to carry more and 
stronger spectral lines than C-H bonds. Figure 4 shows 
the amplitudc-FFT of the length of one C-C bond in the 
LM molecule, more than 220 ps after excitation. The time 
length of the data sample was 65 ps, so that the frequency 
resolution is 1 cm -1 (with Hanning-filtering). Most lines 
are indeed about this wide, and roughly retain their (sim- 
ple) shape but vary in strength to some extent because 
they exchange energy with distant modes. Some are so 
close that they partly blend together. A case in point is the 
C-H stretch band, which extends from 2920 to 2970 cm -1 
(somewhat below experimental value, because of imperfec- 
tions in the model force fields). Since there are about 250 
H atoms in this molecule, we expect about that number 
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Fig. 4. The frequency spectrum of the fluctuations of the 
length of one of the C-C bonds, more than 221 ps after 
excitation. Note the narrowness of the lines, limited by 
the sample length (65 ps). 



of C-H str modes, all in that range, but slightly displaced 
from one another because of differences of environments. 
Clearly, these modes are close enough to strongly interact 
and share energy, like nearly resonant coupled oscillators. 
This is indeed what transpires from the comparison of 
successive FFT's over, say, a 16 ps interval each: the band 
profile changes but the integrated "intensity" remains es- 
sentially constant. This is an example of selectivity in IVR 
(Felker and Zewail 28 ). In this molecule, however, we found 
no case of strong resonance entailing "thermalization" be- 
tween the corresponding oscillators, as shown by Ford et 
al. 16 . By contrast, we did simulate such a resonance by 
exciting some of the stretches of the highly symmetrical 
benzene molecule, and observing the subsequent fluctu- 
ating energy distribution among the six C-H bonds (the 
resonant oscillators) in steady state. 

Overall, spectra like the one in fig. 4 are not grassy 
but consist of distinct, narrow bands whose positions do 
not vary notably in time; the Fourier transforms of the 
corresponding autocorrelation functions deliver identical 
spectra, frequencywise. The motion can therefore be con- 
sidered regular overall. This conclusion is confirmed by 
the autocorrelation functions: e.g., in fig. 5, for the same 
C-C bond length as in fig. 4; note the very strong and 
persistent correlations; the superimposed modulation in- 
dicates fluctuations of line intensities due to energy ex- 
changes between oscillators. Similar results are obtained 
for C-H bonds. Thus, despite heavy spectral congestion, 
the molecule as a whole does not qualify as a thermal- 
izing bath to any individual mode. Reasons for this are 
probably the amorphicity of the structure, which hinders 
close resonances, and the lack of a "phonon continuum" at 
the lowest frequencies of the spectrum (see fig. 6 of 24 ) as 
opposed to the case of molecular crystals 29 , for instance. 
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Fig. 5. The autocorrelation function of the fluctuations in 
length of the same C-C bond as in fig. 4, and during the 
same time. 

For our present purposes, it is important to note that 
this steady state vibrational energy distribution is not 
very sensitive to the details of the H atom capture or, 
in terms of the procedure adopted here, to the details of 
the initial C-H bond perturbation which is meant to sim- 
ulate that capture, nor to the particular C-H bond that is 
perturbed, as long as the excess energy thus deposited re- 
mains of order 4.5 eV. From a more fundamental point of 
view, however, it is interesting to observe that, because of 
the non-linearity of the system, the steady state depends 
on that excess energy as well as on the type of bond which 
is perturbed (C-C, O-H, stretch, bend, torsion, etc.). None 
the less, even if the perturbation consists in altering the 
torsion angle of a C-O-H group, thus exciting one of the 
lowest-frequency modes, vibrational energy finds its way 
up the frequency ladder and most of the normal modes 
are ultimately activated. Figure 6 illustrates the relaxation 
of LM to equilibrium after the torsion angles of six such 
groups were altered by about 100 degrees. Note the dif- 
ferences between the steady state spectrum reached here 
(curve 6) and the one that was reached after altering the 
length of a C-H bond (curve 7, identical with curve 7 of 
fig. 3): the higher frequencies are distinctly weaker here, 
and the opposite is true below 2000 cm -1 . The same trend 
is observed, but to a lesser degree, if energy is instead ini- 
tially deposited in a C-C bond. This sensitivity of the final 
state to the initial perturbation is an indication that sta- 
tistical (thermodynamic) relaxation did not take place. 

It can be shown 24 that the expected spectral density 
of IR radiation from the excited molecule is proportional 
to the product of the square of the frequency with the 
power spectral density of the dipole moment. The lat- 
ter is the FFT of the sum of the autocorrelation func- 
tions of the fluctuations of the dipole moment compo- 
nents about their respective means. Finally, these fluc- 
tuations can be recorded during a run of the code well 
after steady state has settled. This procedure was applied 
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Fig. 6. The relaxation of the cumulative spectral distri- 
bution, Y^kin{v), after perturbation of the torsion angles 
of several C-O-H groups, which adds about 0.6 eV to the 
ground state of the molecule. Spectra 1 to 6 were recorded 
after 8, 16, 65, 200, 450 and 500 ps, respectively. Spectrum 
7 is identical with 7 of fig. 3 and is added here for com- 
parison. Note the general evolution from left to right by 
contrast with fig. 3; also, the asymptotic spectrum 6 is 
poorer than 7 in higher frequencies. 

to the steady state corresponding to fig. 3 (curve 7) and 
fig. 4, and the result is displayed in fig. 7. This shows that 
only a fraction of the vibrational modes is "infrared ac- 
tive" , but one can still recognize the C-H stretch bands 
(aromatic and aliphatic, around 3000 cm -1 ), the in-plane 
C-H bends (~ 1000-1200 cm" 1 ) and the C-H out-of-plane 
bends (~ 750-1000 cm -1 ), and far infrared bands mainly 
associated with oxygen atoms in various environments 23 . 
It should be stressed again that the present computations 
used the MM+ molecular dynamics software, which is 
known to yield only approximate frequencies and intensi- 
ties. A more accurate spectrum can be obtained using the 
PM3 semi-empirical software, albeit on a smaller molecule 
(100 atoms), but still representative of the present amor- 
phous structure 24 . 

5. Chaotic motions 

The conclusions drawn above are strengthened by com- 
parison with conspicuously chaotic motions. Such motions 
are expected when the molecule is immersed in a thermal 
bath. The Hyperchem package allows the simulation of 
this situation by means of Langevin dynamics, in which 
the usual equations of motion are supplemented by terms 
representing both the friction against the bath molecules 
and the random changes in position and velocity upon col- 
lisions with these molecules (Allen and Tildesley 30 ). The 
user may specify the temperature, T, and the friction co- 
efficient, 7; the stronger the latter the tighter the coupling 
of the molecule to the bath. For purposes of comparison 
with fig. 5, we set 7=100 ps -1 and T=35 K, a temperature 



Fig. 7. The infrared spectrum of the isolated molecule LM 
in steady state after chemical excitation by capture of a 
H atom. 
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Fig. 8. Autocorrelations of time fluctuations of the C-C 
bond length for chaotic motions: a) molecule in a ther- 
mal bath at 35 K, simulated by Langevin dynamics (curve 
shifted by (+l,+0.5) for clarity); b) isolated but highly ex- 
cited molecule (12.5 eV). The autocorrelations are char- 
acteristically weaker than in fig. 5. 



at which the energy content of the molecule is about the 
same as above. Figure 8a shows, under these conditions, 
the autocorrelation function of the length of the same C- 
C bond as in fig. 5. The difference is striking: very weak 
correlation, with no regular pattern. 

Theoretical studies also show that chaos is bound to 
set in beyond a total energy threshold, even in an iso- 
lated system. To explore this route to chaos, we first ran 
Langevin dynamics at 3000 K for 0.3 ps. This dumped 
an energy of 12.5 eV uniformly in the molecule. The latter 
was then isolated from the bath and molecular dynamics 
was run in vacuum for 16.5 ps. The corresponding cor- 
relation function is shown as curve b in fig. 8, again for 
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Fig. 9. The cumulative spectral energy distributions a,b 
and c, corresponding, respectively, to autocorrelation 
curve a of fig. 8 (Langevin dynamics), curve b of fig. 8 
(violently excited isolated molecule), and curve 7 of fig. 
3 (mildly excited isolated molecule) . Note the progression 
from quasi-regular to fully chaotic motion, going from c 
to a. 

the same C-C bond length as above. The weak correla- 
tions and highly disordered modulation indicate the par- 
tial onset of chaos, but this cannot be carried much further 
because the molecules dissociate at higher excitation en- 
ergies. Chaos is much easier to induce that way in a small 
and/or symmetric structure like benzene; I confirmed this 
using the same code. 

The frequency spectra of these chaotic motions are also 
distinctly different from the regular spectra: they are very 
grassy and individual modes are no longer recognizable be- 
cause they blend together, even at the highest resolution. 
In the extreme case of Langevin dynamics, all features are 
drowned in the grass which extends to much higher fre- 
quencies. This is apparent from the cumulative spectral 
kinetic energy distributions in fig. 9, where curves 9a, b 
correspond, respectively, to 8a,b while 9c is taken from 
fig. 3, the steady state regular motion. In such graphs, of 
course, the grass is not apparent because it is integrated 
out, so that the progression to chaos from c to a is made 
clearer. 

6. Conclusion 

Given the size and structure of our system, and the type 
and strength of the activation process, there is no theoret- 
ical way to determine, beforehand, the nature of the sub- 
sequent relaxation. The computations performed in this 
study help drawing a realistic classical model of the dy- 
namics to be expected. 

Many markers point to a regular (quasi-periodic) mo- 
tion in steady state: the steady state depends on the initial 
perturbation (although not wildly); the spectral lines are 
narrow; correlations are strong and long lived; the spec- 



tral energy distribution (energy spectrum) is very different 
than for a molecule in a (model) thermal bath; in particu- 
lar, there is no continuum here. However, this steady state 
is not perfectly coherent: the spectral lines vary erratically 
in intensity (although their short term averages are con- 
stant), i.e. there are random, reversible, energy exchanges 
between modes; and the combination frequencies are ab- 
sent from the energy spectrum, i.e. there are no apparent 
phase relations between modes. 

These results are all understandable in terms of more 
or less strongly and non-linearly coupled anharmonic os- 
cillators (see 11,16 and references therein). They fit neatly 
in the picture of partial chaos or broken ergodicity 17 : the 
phase point can wander randomly, but only in restricted 
regions of the energetically available phase space. Such sit- 
uations occur when the vibrational spectrum is too sparse 
and/or the energy content is too low for multiple overlaps 
to occur between non-linear resonances. In the present in- 
stance, this is likely due to the limited size, and to the 
amorphicity and inhomogeneity of the structure. 

From comparisons between corresponding classical and 
quantal motions (see Noid et al. 19 ), it appears that quan- 
tum chaos does not set in before classical chaos is well es- 
tablished. It seems, therefore, that, even under the more 
realistic, but "corresponding", laws of quantum mechan- 
ics, the vibrational motions of our activated molecule will 
not reach statistical (thermodynamic) equilibrium. If that 
is the case, then the classical spectral intensity (fig. 6) 
may be interpreted as the ensemble average quantal en- 
ergy content at the same frequency. One can therefore ex- 
pect that such molecules, isolated and chemically excited, 
will ultimately deliver a significant part of their initial en- 
ergy content in the form of near and mid-IR radiation. 

The question of the upper size limit beyond which 
these conclusions no longer apply remains open to dis- 
cussion. It is known that the vibrational energy de- 
posited in a molecular crystals will ultimately flow into the 
phonon bath 29 ' 31 . Phonons are very low-frequency (<I00 
cm -1 ) plane propagating waves which can transport en- 
ergy through the structure. They can interact upon hitting 
a structural or chemical defect, and exchange energy like 
perfect gas particles in a vessel, so that thermal equilib- 
rium can set in. But, for the localized vibrations (vibrons), 
which we have been considering above, to be efficiently 
converted into phonons, the spectrum of the latter must 
be dense enough as to form a quasi-continuum to which 
the Fermi Golden Rule of transitions can be applied. This 
is only possible if the structure is regular and infinite, 
like ideal crystals; very irregular structures obviously can 
hardly support phonons, a fact which translates into very 
low thermal conductivities and long thermal relaxation 
times 32 . Indeed, the relaxation times of some glasses at 
low temperatures have been found to range in the hours 33 . 
It is likely, therefore, that a grain of our model material 
will have to be quite large before it can give rise to a 
phonon spectrum of sufficient density that it constitutes 
an efficient thermal bath. 
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